#ifndef LEASTSQUARESOLVER_HPP
#define LEASTSQUARESOLVER_HPP

#include "matrix.hpp"
#include "matrixdecomposition.hpp"
#include <iostream>
namespace TermStructure 
{
	//! Least Square Solver based on two matrices.
	/*! Solve for lease square fit, the solver is LUDEcomposition Solver
	*/
	template <typename T, typename TSize = long>
	class LeastSquareSolver
	{
	public:
		//! Fit the Square
		/*! \param input the input matrix.
		    \param load the load.
			\param result where the result will be stored. Please pre-locate the space
		*/
		template<typename TInputMatrix, typename TLoadMatrix, typename TResultMatrix>
		static void Fit(const TInputMatrix &input, const TLoadMatrix &Load, TResultMatrix &result)
		{
			SimpleMatrix<T,TSize> left(input.GetLength(1),input.GetLength(1)), temp(input.GetLength(1),input.GetLength(1));
			SimpleMatrix<T,TSize> right(input.GetLength(1));
			for(TSize i = 0; i < input.GetLength(1); i++)
			{
				for(TSize j = 0; j <= i; j++)
				{
					left(i,j) = 0;
					for(TSize k = 0; k < input.GetLength(0); k++)
					{
						left(i,j) += input(k,i)*input(k,j);
					}
					left(j,i) = left(i,j);
				}
			}
			for(TSize i = 0; i < input.GetLength(1); i++)
			{
				right(i,0) = 0;
				for(TSize k = 0; k < input.GetLength(0); k++)
				{
					right(i,0) += input(k,i)*Load(k,0);
				}
			}
			LUDecompositionSolve(left, right, result, temp);
		}
	};
}


#endif
